Hands-on_Ex10

Introduction

This project aims to explore various spatial interaction models (SIMs) using R, following the structure provided in the linked materials. We will process and visualize flow data using an Origin-Destination (OD) data set, build an OD matrix, and analyze passenger volume by origin and destination bus stops, following the guidelines from Chapter 15.

Libraries and Setup

First, we need to load the necessary libraries for our analysis.

# Loading required libraries
pacman::p_load(tmap, sf, DT, stplanr, tidyverse)

Preparing the Flow Data

Importing the OD Data

We start by importing the Passenger Volume by Origin Destination Bus Stops data set downloaded from LTA DataMall.

# Importing OD data
odbus <- read_csv("/Users/sharon/OneDrive - Singapore Management University/ISSS626 Data/Hands-on_Ex/Hands-on_Ex10/data/aspatial/origin_destination_bus_202210.csv")

Rows: 5122925 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (3): YEAR_MONTH, DAY_TYPE, PT_TYPE
dbl (4): TIME_PER_HOUR, ORIGIN_PT_CODE, DESTINATION_PT_CODE, TOTAL_TRIPS

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# Displaying the structure of OD data
glimpse(odbus)

Rows: 5,122,925
Columns: 7
$ YEAR_MONTH          <chr> "2022-10", "2022-10", "2022-10", "2022-10", "2022-…
$ DAY_TYPE            <chr> "WEEKDAY", "WEEKENDS/HOLIDAY", "WEEKENDS/HOLIDAY",…
$ TIME_PER_HOUR       <dbl> 10, 10, 7, 11, 16, 16, 20, 7, 7, 11, 11, 8, 11, 11…
$ PT_TYPE             <chr> "BUS", "BUS", "BUS", "BUS", "BUS", "BUS", "BUS", "…
$ ORIGIN_PT_CODE      <dbl> 65239, 65239, 23519, 52509, 54349, 54349, 43371, 8…
$ DESTINATION_PT_CODE <dbl> 65159, 65159, 23311, 42041, 53241, 53241, 14139, 9…
$ TOTAL_TRIPS         <dbl> 2, 1, 2, 1, 1, 4, 1, 3, 1, 5, 2, 5, 15, 40, 1, 1, …

# Convert numeric to character data type
odbus$ORIGIN_PT_CODE <- as.factor(odbus$ORIGIN_PT_CODE)
odbus$DESTINATION_PT_CODE <- as.factor(odbus$DESTINATION_PT_CODE)

Extracting Study Data

To extract commuting flows during weekday mornings (6-9 am):

odbus6_9 <- odbus %>%
  filter(DAY_TYPE == "WEEKDAY") %>%
  filter(TIME_PER_HOUR >= 6 & TIME_PER_HOUR <= 9) %>%
  group_by(ORIGIN_PT_CODE, DESTINATION_PT_CODE) %>%
  summarise(TRIPS = sum(TOTAL_TRIPS))

`summarise()` has grouped output by 'ORIGIN_PT_CODE'. You can override using
the `.groups` argument.

# Save the extracted data for future use
write_rds(odbus6_9, "/Users/sharon/OneDrive - Singapore Management University/ISSS626 Data/Hands-on_Ex/Hands-on_Ex10/data/rds/odbus6_9.rds")

# Display the data in a table format
datatable(odbus6_9)

Warning in instance$preRenderHook(instance): It seems your data is too big for
client-side DataTables. You may consider server-side processing:
https://rstudio.github.io/DT/server.html

Working with Geospatial Data

Importing Geospatial Data

We will use two geospatial data sets: BusStop and MPSZ-2019, both in ESRI shapefile format.

# Import BusStop shapefile
busstop <- st_read(dsn = "/Users/sharon/OneDrive - Singapore Management University/ISSS626 Data/Hands-on_Ex/Hands-on_Ex10/data/geospatial", layer = "BusStop") %>%
  st_transform(crs = 3414)

Reading layer `BusStop' from data source 
  `/Users/sharon/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/ISSS626 Data/Hands-on_Ex/Hands-on_Ex10/data/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 5159 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 3970.122 ymin: 26482.1 xmax: 48284.56 ymax: 52983.82
Projected CRS: SVY21

# Import MPSZ-2019 shapefile
mpsz <- st_read(dsn = "/Users/sharon/OneDrive - Singapore Management University/ISSS626 Data/Hands-on_Ex/Hands-on_Ex10/data/geospatial", layer = "MPSZ-2019") %>%
  st_transform(crs = 3414)

Reading layer `MPSZ-2019' from data source 
  `/Users/sharon/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/ISSS626 Data/Hands-on_Ex/Hands-on_Ex10/data/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 332 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS:  WGS 84

# Save MPSZ data for future use
write_rds(mpsz, "/Users/sharon/OneDrive - Singapore Management University/ISSS626 Data/Hands-on_Ex/Hands-on_Ex10/data/rds/mpsz.rds")

Geospatial Data Wrangling

We will populate the planning subzone code of mpsz into busstop and merge it with the OD data.

# Intersection between busstop and mpsz
busstop_mpsz <- st_intersection(busstop, mpsz) %>%
  select(BUS_STOP_N, SUBZONE_C) %>%
  st_drop_geometry()

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

# Save the resultant data
write_rds(busstop_mpsz, "/Users/sharon/OneDrive - Singapore Management University/ISSS626 Data/Hands-on_Ex/Hands-on_Ex10/data/rds/busstop_mpsz.rds")

# Append subzone codes to OD data
od_data <- left_join(odbus6_9, busstop_mpsz, by = c("ORIGIN_PT_CODE" = "BUS_STOP_N")) %>%
  rename(ORIGIN_BS = ORIGIN_PT_CODE, ORIGIN_SZ = SUBZONE_C, DESTIN_BS = DESTINATION_PT_CODE)

Warning in left_join(odbus6_9, busstop_mpsz, by = c(ORIGIN_PT_CODE = "BUS_STOP_N")): Detected an unexpected many-to-many relationship between `x` and `y`.
ℹ Row 55491 of `x` matches multiple rows in `y`.
ℹ Row 161 of `y` matches multiple rows in `x`.
ℹ If a many-to-many relationship is expected, set `relationship =
  "many-to-many"` to silence this warning.

# Handle duplicates
od_data <- unique(od_data)
od_data <- left_join(od_data, busstop_mpsz, by = c("DESTIN_BS" = "BUS_STOP_N"))

Warning in left_join(od_data, busstop_mpsz, by = c(DESTIN_BS = "BUS_STOP_N")): Detected an unexpected many-to-many relationship between `x` and `y`.
ℹ Row 74 of `x` matches multiple rows in `y`.
ℹ Row 1379 of `y` matches multiple rows in `x`.
ℹ If a many-to-many relationship is expected, set `relationship =
  "many-to-many"` to silence this warning.

od_data <- unique(od_data)
od_data <- od_data %>%
  rename(DESTIN_SZ = SUBZONE_C) %>%
  drop_na() %>%
  group_by(ORIGIN_SZ, DESTIN_SZ) %>%
  summarise(MORNING_PEAK = sum(TRIPS))

`summarise()` has grouped output by 'ORIGIN_SZ'. You can override using the
`.groups` argument.

# Save the final OD data
write_rds(od_data, "/Users/sharon/OneDrive - Singapore Management University/ISSS626 Data/Hands-on_Ex/Hands-on_Ex10/data/rds/od_data_fii.rds")

Visualising Spatial Interaction

Removing Intra-Zonal Flows

We will exclude intra-zonal flows before creating desire lines.

od_data_fij <- od_data[od_data$ORIGIN_SZ != od_data$DESTIN_SZ, ]

# Save the filtered data
write_rds(od_data_fij, "/Users/sharon/OneDrive - Singapore Management University/ISSS626 Data/Hands-on_Ex/Hands-on_Ex10/data/rds/od_data_fij.rds")

Creating Desire Lines

Next, we create the desire lines using the stplanr package.

# Creating desire lines
flowLine <- od2line(flow = od_data_fij, zones = mpsz, zone_code = "SUBZONE_C")

Creating centroids representing desire line start and end points.

# Save desire lines for future use
write_rds(flowLine, "/Users/sharon/OneDrive - Singapore Management University/ISSS626 Data/Hands-on_Ex/Hands-on_Ex10/data/rds/flowLine.rds")

Visualising the Desire Lines

We visualize the resulting desire lines using tmap.

# Visualising all desire lines
tm_shape(mpsz) +
  tm_polygons() +
  tm_shape(flowLine) +
  tm_lines(lwd = "MORNING_PEAK", style = "quantile", scale = c(0.1, 1, 3, 5, 7, 10), n = 6, alpha = 0.3)

Warning in g$scale * (x/maxW): longer object length is not a multiple of
shorter object length

# Focus on flows greater than 5000
tm_shape(mpsz) +
  tm_polygons() +
  tm_shape(flowLine %>% filter(MORNING_PEAK >= 5000)) +
  tm_lines(lwd = "MORNING_PEAK", style = "quantile", scale = c(0.1, 1, 3, 5, 7, 10), n = 6, alpha = 0.3)

Warning in g$scale * (x/maxW): longer object length is not a multiple of
shorter object length

Gravity Model Calibration

To analyze and predict interactions between origins and destinations, we apply a gravity model calibration. This model helps understand the influence of spatial attributes on trip volume.

library(MASS)

Attaching package: 'MASS'

The following object is masked from 'package:dplyr':

    select

# Define a gravity model using MORNING_PEAK as the dependent variable
gravity_model <- glm(MORNING_PEAK ~ ORIGIN_SZ + DESTIN_SZ, family = poisson, data = od_data_fij)


# Check if the model converged
if (!is.null(gravity_model)) {
  # Summary of the gravity model
  summary(gravity_model)
  # Save the model for future use
  tryCatch({
    saveRDS(gravity_model, "/Users/sharon/OneDrive - Singapore Management University/ISSS626 Data/Hands-on_Ex/Hands-on_Ex10/data/rds/gravity_model.rds")
  }, error = function(e) {
    message("Error saving model: ", e$message)
  })
} else {
  message("Gravity model did not converge. Please check the data or try a different model specification.")
}